#include<iostream>
#include<vector>
#include<queue>
#include<malloc.h>
#include <math.h>

void A_B(int order)
  
{
  int i,n,q;
  double u,k1,k2,T1,T2,y0,y1,y2,y3,z0,z1,z2,z3,b1,b2,b3,b0;
  double *u1,*u2,*u4,*u5;
  double e1,e2,emaxnorm,premaxnorm,Crate;
  double start,finish,duration;

  u1=(double *)malloc(sizeof(double)*10000000);
  u2=(double *)malloc(sizeof(double)*10000000);
  u4=(double *)malloc(sizeof(double)*10000000);
  u5=(double *)malloc(sizeof(double)*10000000);
  u1[0]=0.994;u2[0]=0;u4[0]=0;u5[0]=-2.0015851063790825224;u=1/81.45;
  T1=17.06521656015796;T2=19.14045706162071;

  if (order==1){
    b0=1.0;b1=0.0;b2=0.0;b3=0.0;}
  else if (order==2){
    b0=-0.5;b1=1.5;b2=0.0;b3=0.0;}
  else if (order==3){
    b0=5.0/12;b1=-16.0/12;b2=23.0/12;b3=0.0;}
  else{
    b0=-9.0/24;b1=37.0/24;b2=-59.0/24;b3=55.0/24;}

   std::cout<<"p="<<order<<std::endl;
std::cout<<"the A_B method(situation1):"<<std::endl;

      std::cout<<"  节点数  "<<"  运行时间  "<<" 误差的最大范数  "<<"  收敛阶  "<<std::endl;

      n=4e3;;k1=T1/n;
        for (i=1;i<=3;i++){
    u1[i]=u1[i-1]+k1*u4[i-1];
    u2[i]=u2[i-1]+k1*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    u4[i]=u4[i-1]+k1*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    u5[i]=u5[i-1]+k1*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
  }
  
  for (i=4;i<=n;i++){ 
    u1[i]=u1[i-1]+k1*(b0*u4[i-1]+b1*u4[i-2]+b2*u4[i-3]+b3*u4[i-4]);
    u2[i]=u2[i-1]+k1*(b0*u5[i-1]+b1*u5[i-2]+b2*u5[i-3]+b3*u5[i-4]);
    y0=u2[i-4]*u2[i-4]+(u1[i-4]+u-1)*(u1[i-4]+u-1);y0=pow(y0,1.5);
    y1=u2[i-3]*u2[i-3]+(u1[i-3]+u-1)*(u1[i-3]+u-1);y1=pow(y1,1.5);
    y2=u2[i-2]*u2[i-2]+(u1[i-2]+u-1)*(u1[i-2]+u-1);y2=pow(y2,1.5);
    y3=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y3=pow(y3,1.5);
    z0=u2[i-4]*u2[i-4]+(u1[i-4]+u)*(u1[i-4]+u);z0=pow(z0,1.5);
    z1=u2[i-3]*u2[i-3]+(u1[i-3]+u)*(u1[i-3]+u);z1=pow(z1,1.5);
    z2=u2[i-2]*u2[i-2]+(u1[i-2]+u)*(u1[i-2]+u);z2=pow(z2,1.5);
    z3=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);z3=pow(z3,1.5);
    u4[i]=u4[i-1]+k1*(b3*(2*u5[i-4]+u1[i-4]-(u*(u1[i-4]+u-1))/y0-((1-u)*(u1[i-4]+u))/z0)+b2*(2*u5[i-3]+u1[i-3]-(u*(u1[i-3]+u-1))/y1-((1-u)*(u1[i-3]+u))/z1)+
		      b1*(2*u5[i-2]+u1[i-2]-(u*(u1[i-2]+u-1))/y2-((1-u)*(u1[i-2]+u))/z2)+b0*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y3-((1-u)*(u1[i-1]+u))/z3));
    u5[i]=u5[i-1]+k1*(b3*(u2[i-4]-2*u4[i-4]-(u*u2[i-4])/y0-((1-u)*u2[i-4])/z0)+b2*(u2[i-3]-2*u4[i-3]-(u*u2[i-3])/y1-((1-u)*u2[i-3])/z1)+
		      b1*(u2[i-2]-2*u4[i-2]-(u*u2[i-2])/y2-((1-u)*u2[i-2])/z2)+b0*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y3-((1-u)*u2[i-1])/z3));
  }

  e1=u1[n]-u1[0];
  e2=u2[n]-u2[0];
  e1=fabs(e1);e2=fabs(e2);
  emaxnorm=std::max(e1,e2);



      //
  for (q=1;q<=10;q++){
  premaxnorm=emaxnorm;
  start = clock();    
  n=n*2;
  k1=T1/n;
  
  for (i=1;i<=3;i++){
    u1[i]=u1[i-1]+k1*u4[i-1];
    u2[i]=u2[i-1]+k1*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    u4[i]=u4[i-1]+k1*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    u5[i]=u5[i-1]+k1*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
  }
  
  for (i=4;i<=n;i++){ 
    u1[i]=u1[i-1]+k1*(b0*u4[i-1]+b1*u4[i-2]+b2*u4[i-3]+b3*u4[i-4]);
    u2[i]=u2[i-1]+k1*(b0*u5[i-1]+b1*u5[i-2]+b2*u5[i-3]+b3*u5[i-4]);
    y0=u2[i-4]*u2[i-4]+(u1[i-4]+u-1)*(u1[i-4]+u-1);y0=pow(y0,1.5);
    y1=u2[i-3]*u2[i-3]+(u1[i-3]+u-1)*(u1[i-3]+u-1);y1=pow(y1,1.5);
    y2=u2[i-2]*u2[i-2]+(u1[i-2]+u-1)*(u1[i-2]+u-1);y2=pow(y2,1.5);
    y3=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y3=pow(y3,1.5);
    z0=u2[i-4]*u2[i-4]+(u1[i-4]+u)*(u1[i-4]+u);z0=pow(z0,1.5);
    z1=u2[i-3]*u2[i-3]+(u1[i-3]+u)*(u1[i-3]+u);z1=pow(z1,1.5);
    z2=u2[i-2]*u2[i-2]+(u1[i-2]+u)*(u1[i-2]+u);z2=pow(z2,1.5);
    z3=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);z3=pow(z3,1.5);
    u4[i]=u4[i-1]+k1*(b3*(2*u5[i-4]+u1[i-4]-(u*(u1[i-4]+u-1))/y0-((1-u)*(u1[i-4]+u))/z0)+b2*(2*u5[i-3]+u1[i-3]-(u*(u1[i-3]+u-1))/y1-((1-u)*(u1[i-3]+u))/z1)+
		      b1*(2*u5[i-2]+u1[i-2]-(u*(u1[i-2]+u-1))/y2-((1-u)*(u1[i-2]+u))/z2)+b0*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y3-((1-u)*(u1[i-1]+u))/z3));
    u5[i]=u5[i-1]+k1*(b3*(u2[i-4]-2*u4[i-4]-(u*u2[i-4])/y0-((1-u)*u2[i-4])/z0)+b2*(u2[i-3]-2*u4[i-3]-(u*u2[i-3])/y1-((1-u)*u2[i-3])/z1)+
		      b1*(u2[i-2]-2*u4[i-2]-(u*u2[i-2])/y2-((1-u)*u2[i-2])/z2)+b0*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y3-((1-u)*u2[i-1])/z3));
  }
  e1=u1[n]-u1[0];
  e2=u2[n]-u2[0];
  e1=fabs(e1);e2=fabs(e2);
  emaxnorm=std::max(e1,e2);
Crate=premaxnorm/emaxnorm;Crate=log2(Crate);
  finish = clock();
  duration = (double)(finish - start) / CLOCKS_PER_SEC;
  std::cout<<"  "<<n<<"     " <<duration<<"    "  <<emaxnorm<<"      "<<Crate  <<std::endl;
 
  }


  //


   T1=T2,u1[0]=0.87978;u2[0]=0;u4[0]=0;u5[0]=-0.3797;
   std::cout<<"the A_B method(situation2):"<<std::endl;
      std::cout<<"  节点数  "<<"  运行时间  "<<" 误差的最大范数  "<<"  收敛阶  "<<std::endl;


      n=5e3;;k1=T1/n;
        for (i=1;i<=3;i++){
    u1[i]=u1[i-1]+k1*u4[i-1];
    u2[i]=u2[i-1]+k1*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    u4[i]=u4[i-1]+k1*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    u5[i]=u5[i-1]+k1*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
  }
  
  for (i=4;i<=n;i++){ 
    u1[i]=u1[i-1]+k1*(b0*u4[i-1]+b1*u4[i-2]+b2*u4[i-3]+b3*u4[i-4]);
    u2[i]=u2[i-1]+k1*(b0*u5[i-1]+b1*u5[i-2]+b2*u5[i-3]+b3*u5[i-4]);
    y0=u2[i-4]*u2[i-4]+(u1[i-4]+u-1)*(u1[i-4]+u-1);y0=pow(y0,1.5);
    y1=u2[i-3]*u2[i-3]+(u1[i-3]+u-1)*(u1[i-3]+u-1);y1=pow(y1,1.5);
    y2=u2[i-2]*u2[i-2]+(u1[i-2]+u-1)*(u1[i-2]+u-1);y2=pow(y2,1.5);
    y3=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y3=pow(y3,1.5);
    z0=u2[i-4]*u2[i-4]+(u1[i-4]+u)*(u1[i-4]+u);z0=pow(z0,1.5);
    z1=u2[i-3]*u2[i-3]+(u1[i-3]+u)*(u1[i-3]+u);z1=pow(z1,1.5);
    z2=u2[i-2]*u2[i-2]+(u1[i-2]+u)*(u1[i-2]+u);z2=pow(z2,1.5);
    z3=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);z3=pow(z3,1.5);
    u4[i]=u4[i-1]+k1*(b3*(2*u5[i-4]+u1[i-4]-(u*(u1[i-4]+u-1))/y0-((1-u)*(u1[i-4]+u))/z0)+b2*(2*u5[i-3]+u1[i-3]-(u*(u1[i-3]+u-1))/y1-((1-u)*(u1[i-3]+u))/z1)+
		      b1*(2*u5[i-2]+u1[i-2]-(u*(u1[i-2]+u-1))/y2-((1-u)*(u1[i-2]+u))/z2)+b0*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y3-((1-u)*(u1[i-1]+u))/z3));
    u5[i]=u5[i-1]+k1*(b3*(u2[i-4]-2*u4[i-4]-(u*u2[i-4])/y0-((1-u)*u2[i-4])/z0)+b2*(u2[i-3]-2*u4[i-3]-(u*u2[i-3])/y1-((1-u)*u2[i-3])/z1)+
		      b1*(u2[i-2]-2*u4[i-2]-(u*u2[i-2])/y2-((1-u)*u2[i-2])/z2)+b0*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y3-((1-u)*u2[i-1])/z3));
  }

  e1=u1[n]-u1[0];
  e2=u2[n]-u2[0];
  e1=fabs(e1);e2=fabs(e2);
  emaxnorm=std::max(e1,e2);



      //
  for (q=1;q<=10;q++){
  premaxnorm=emaxnorm;
  start = clock();    
  n=n*2;
  k1=T1/n;
  
  for (i=1;i<=3;i++){
    u1[i]=u1[i-1]+k1*u4[i-1];
    u2[i]=u2[i-1]+k1*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    u4[i]=u4[i-1]+k1*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    u5[i]=u5[i-1]+k1*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
  }
  
  for (i=4;i<=n;i++){ 
    u1[i]=u1[i-1]+k1*(b0*u4[i-1]+b1*u4[i-2]+b2*u4[i-3]+b3*u4[i-4]);
    u2[i]=u2[i-1]+k1*(b0*u5[i-1]+b1*u5[i-2]+b2*u5[i-3]+b3*u5[i-4]);
    y0=u2[i-4]*u2[i-4]+(u1[i-4]+u-1)*(u1[i-4]+u-1);y0=pow(y0,1.5);
    y1=u2[i-3]*u2[i-3]+(u1[i-3]+u-1)*(u1[i-3]+u-1);y1=pow(y1,1.5);
    y2=u2[i-2]*u2[i-2]+(u1[i-2]+u-1)*(u1[i-2]+u-1);y2=pow(y2,1.5);
    y3=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y3=pow(y3,1.5);
    z0=u2[i-4]*u2[i-4]+(u1[i-4]+u)*(u1[i-4]+u);z0=pow(z0,1.5);
    z1=u2[i-3]*u2[i-3]+(u1[i-3]+u)*(u1[i-3]+u);z1=pow(z1,1.5);
    z2=u2[i-2]*u2[i-2]+(u1[i-2]+u)*(u1[i-2]+u);z2=pow(z2,1.5);
    z3=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);z3=pow(z3,1.5);
    u4[i]=u4[i-1]+k1*(b3*(2*u5[i-4]+u1[i-4]-(u*(u1[i-4]+u-1))/y0-((1-u)*(u1[i-4]+u))/z0)+b2*(2*u5[i-3]+u1[i-3]-(u*(u1[i-3]+u-1))/y1-((1-u)*(u1[i-3]+u))/z1)+
		      b1*(2*u5[i-2]+u1[i-2]-(u*(u1[i-2]+u-1))/y2-((1-u)*(u1[i-2]+u))/z2)+b0*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y3-((1-u)*(u1[i-1]+u))/z3));
    u5[i]=u5[i-1]+k1*(b3*(u2[i-4]-2*u4[i-4]-(u*u2[i-4])/y0-((1-u)*u2[i-4])/z0)+b2*(u2[i-3]-2*u4[i-3]-(u*u2[i-3])/y1-((1-u)*u2[i-3])/z1)+
		      b1*(u2[i-2]-2*u4[i-2]-(u*u2[i-2])/y2-((1-u)*u2[i-2])/z2)+b0*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y3-((1-u)*u2[i-1])/z3));
  }
  e1=u1[n]-u1[0];
  e2=u2[n]-u2[0];
  e1=fabs(e1);e2=fabs(e2);
  emaxnorm=std::max(e1,e2);
Crate=premaxnorm/emaxnorm;Crate=log2(Crate);
  finish = clock();
  duration = (double)(finish - start) / CLOCKS_PER_SEC;
  std::cout<<"  "<<n<<"     " <<duration<<"    "  <<emaxnorm<<"      "<<Crate  <<std::endl;
 
  }

      
}
